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Abstract. We study a Bose-Fermi mixture within the framework of the mean- 
field theory, including three possible regimes for the fermionic species: fully 
polarized, BCS, and unitarity. Starting from the 3D description and using the 
variational approximation (VA), we derive ID and 2D systems of equations, under 
the corresponding confining potentials. This method produces a pair of nonlinear 
Schrodinger (NTS) equations coupled to algebraic equations for the transverse widths 
of the confined state. The equations incorporate interactions between atoms of the 
same species and between the species, assuming that the latter can be manipulated 
by means of the Feshbach resonance (FR). As an application, we explore spatial 
density correlations in the ground state (GS) between the species, concluding that 
they strongly depend on the sign and strength of the inter-species interaction. Also 
studied are the dynamics of the mixture in a vicinity of the GS and the corresponding 
spatiotemporal inter-species correlation. The correlations are strongly affected by the 
fermionic component, featuring the greatest variation in the unitary regime. Results 
produced by the VA are verified by comparison with full numerical solutions. 
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1. Introduction 

Cold atom gases trapped in magneto-optical potentials have tremendously extended 
the range of experimental possibilities in condensed-matter physics. The possibility of 
controlling the intensity of the interactions, the depth of the trapping potential lattices, 
and their geometry enables emulating a broad variety of solid-state settings and the 
creation of novel states of the quantum matter [l]-[8]. 

In this context, a great variety of experimental and theoretical studies were 
stimulated by the possibility of mixing degenerate bosonic and fermionic gases. The 
hrst experimental works on Bose-Fermi mixtures (BFMs) were conducted with lithium 
isotopes [9, 10], which was a precursor to the study of new combinations, such as 
^^■^Yb-^Li [11] and ®'^Rb-‘^°K [12]. Particular attention has been paid to heavy-atom 
mixtures, e.g., ®^Sr-®^Sr [13] isotopes with a large nuclear spin, that have been proposed 
as prototypes for handling the quantum information. One of important points of the 
work in this direction is the use of Feshbach resonances (FRs) in the mixtures, as 
they make it possible to control the two-body interactions between the species. The 
FR has been observed for the ®^Rb-'^*’K mixture [14, 15], and a giant FR has been 
reported for ®^Rb-®Li [16]. Further, hve FRs were observed in ®Li-^^^Cs [17], and over 
30 resonances are known to occur in the ^^Na-^'^K mixture [18]. In Ref. [19], a multiple 
heteronuclear FRs was reported in a triply quantum-degenerate mixture of bosonic 
and two fermionic species, and ®Li. In the latter context, a wide s-wave resonance 
for combination transforming the system into a strongly interacting isotopic 

BFM immersed into a Fermi sea. 

Parallel to the experiment, several theoretical works addressed the dynamics of 
BFMs [20]-[26]. An approach that has proven to be very useful for describing the 
ground state (GS) of the mixtures is provided by the mean held theory [27]-[31]. In this 
context, the FRs for mixtures of ^^Na-®Li, ®^Rb-"^°K, ®^Rb-®Li, ^He-"^He, ^^^Yb-^^'^Yb, 
and ®’^Sr-®'^Sr have been studied [32, 33]. 

The possibility of holding ultracold gases in magneto-optical traps makes it possible 
to create effectively one- or two-dimensional (ID or 2D) situations. A great number of 
works addressed such ID settings [34]-[41]. In particular, the variational approximation 
makes it possible to reduce the 3D description to ID or 2D [42, 43]. In many cases, 
these approximations do not take into regard the fact that the transverse conhnement 
width(s) may vary with the particle density in the unconhned direction. The variational 
approximation that does take this fact into account was developed for bosons in Ref. 
[44] and for fermions in Ref. [45], demonstrating the high accuracy of the ID and 2D 
approximations in comparison with the full the 3D description. 

In this paper we analyze the BFM at zero temperature conhned to ID or 2D 
by external potentials. First, we consider the cigar-shaped conhguration produced by 
the strong conhnement with the cylindrical symmetry acting in the transverse plane, 
while a generic weak potential may act in the axial direction. Second, we consider 
the disc-shaped BFM corresponding to the strong conhnement applied in one direction. 
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and a weak 2D potential acting in the perpendicular plane. We use the variational 
method developed in Refs. [44, 45] to derive effective ID and 2D nonlinear Schrodinger 
(NLS) equations for the cigar-shaped and disc-shaped conhgurations, respectively. In 
particular, our results apply to the ^Li-®Li BFM. In both cases (ID and 2D), we 
analyze the influence of the density of both species on the width of the mixture in 
the conhned dimensions, the main control parameter being the inter-species scattering 
length, obf- The analysis is performed for three different regimes of the fermionic 
component: polarized, BCS, and unitarity. Using the dimensionally-reduced system, we 
consider spatial correlations between densities of particles of both species. Considering 
perturbations around the GS, we study the spatiotemporal correlation between the 
species too. 

The paper is organized as follows: In Sec. 2, we derive ID and 2D variational 
equations from the full 3D description. In Sec. 3, we apply our equations to the numerical 
analysis of the GS in the BFM, and consideration of the spatial correlation between the 
species. For the ID setting, we calculate spatiotemporal correlation in the form of the 
so-called Pearson coefficient, considering perturbations around the GS. Conclusions are 
presented in Sec. 4. 


2. The variational approximation 


We consider a diluted superfluid mixture formed by A^b bosonic atoms of mass 
mB, and iVp fermionic atoms of mass mp and spin sp. The atoms interact 
through the pseudopotential, (5(r) [1]. We assume that bosons form a Bose-Einstein 
condensate (BEG), described by the Gross-Pitaevskii equation [1], while the local 
density approximation [1] applies to the description of the weakly interacting fermionic 
component. Accordingly, the dynamical equations for the BFM can be derived from 
action S, 


S — dtdv (/2b T /2f + /2bf) 


( 1 ) 


where Cb and C-p are the Lagrangian densities of the Bose and Fermi components, while 
/2 bf accounts for the interaction between the them: 


/2b — 


ih 
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|VTb|'-Ub|Tb|'--^?b|Tb|",(2) 
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( 3 ) 


( 4 ) 


Here Gp = [bvr^/ (2 sf + is a constant that depends on spin sp; qb = ^Trhd’aB/rnp,, 
gp = 47rh^(aF/mF)[2S'F/(2S'F -|- 1)], and ^^bf = ^'wh^app/mpp are the three interaction 
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Table 1. Ai, A 2 , and sp for three different regimes of the Fermi gas. 


Regime 

Ai 

A 2 


sf 

Polarized 

1 

1 

1 

0 

BCS 

2 

4 

1 

1/2 

Unitarity 

2 

4 

0.44 

1/2 


parameters of the mixture, with ob, of and obf being the respective scattering lengths; 
tubf = mBm-p/{mB + m-p) is the reduced mass; and U-q/p are external potentials acting 
on bosons/fermions. Complex wave functions ft^n/F (r,t) are normalized to the numbers 
of particles, A^b/f- Parameters Ai, A 2 , and ^ in the fermionic Lagrangian density (3), 
along with sp, determine three different regimes [46]-[49] listed in the Table 1. 

We apply the formalism developed below to the ^Li-®Li mixture, with the same 
scattering parameter for both species, ob/f = 5nm. The use of isotopes of the same 
alkali element is suggested by the similarity of their electric polarizability, thus implying 
similar external potentials induced by an optical trap. Unless specihed otherwise, in 
what follows below we consider conhgurations with fully polarized fermions. Values of 
parameters for all fermionic regimes are collected in Table 1. The BCS and unitarity 
regimes involve more than one spin state of fermions, hence the magnetic trap will split 
the respective spin energy levels. For this reason, we assume the presence of the optical 
trap, which supports equal energy levels for all the spin states, making it possible to 
discriminate different regimes of the interaction in the Bose-Fermi mixture. In the BCS 
and unitarity regimes, we assume balanced populations of the two spin components. 

Finally, varying action S with respect to Tg and to Tp, we derive the following 
system of nonlinear Schrodinger equations for bosons and fermions: 


ihdt^B = 


-|- QbY ^+ 5'BF|fkF|^ + Ub 


2m B 




B, 


( 5 ) 


ih 


dt'^p — 




2X2mp 


V^-|-5'f|^f| +fl'BF|^s| + Up 


2m. 


^C'fI^i 


14/3 




( 6 ) 


As is known from previous works, the mean-held description of the fermionic 
component, based on the effective NFS equation (6) is valid in the hydrodynamic 
approximation, which implies that the Fermi component forms a correlated superhuid. 
In static situations, this approximation is relevant for sufhciently smooth conhgurations, 
whose characteristic spatial scale is much greater than the de Broglie wavelength at the 
Fermi surface [50, 51]. In various dynamical settings, both pure fermionic and mixed 
Bose-Fermi ones, the hydrodynamic approach remains relevant for slow quasi-adiabatic 
evolution [42]-[49]; [52, 53]. 

While numerical integration of this system in the 3D form system is very heavy 
computationally, the ehective dimension may be reduced to ID or 2D when the system 
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is tightly confined by a trapping potential. To this end, the variational method is 
employed, making use of the factorization of the 3D wave function, which includes a 
Gaussian ansatz in the tightly confined transverse directions. The factorization has been 
widely used for Bose and Fermi system separately, as it shown in Refs. [44] and [45], 
respectively. In next two sections we apply the variational technique to a Bose-Fermi 
mixture, obtaining results that agree well with full 3D simulations. On the other hand, 
the applicability of the Gaussian ansatz for approximating of the transverse part of the 
3D wave function is determined by a relation of the transverse-confinement strength, and 
the effective strength of the bosonic or fermionic nonlinearity. In case the latter factor 
is stronger, the transverse wave function may be used in the form of the Thomas-Fermi 
(TF) approximation [53]. 


2.1. The one-dimensional variational system 


The common form of the confinement is provided by the tight harmonic-oscillator 
potential acting in the trapping dimensions. Thus, the confinement to ID (the cigar¬ 
shaped configuration along the z-axis) corresponds to the following potentials, which 
include weak axial components, Vid,b/f: 

Vb/F (r, t) = -"iB/Fl^i,B/F + y^) + ^1D,B/F {z, t) . (7) 

Assuming that the transverse trapping potential is strong enough, the dimensional 
reduction is carried out by means of the usual factorized ansatz for the wave 
functions, 4/B/F(r, t) = vb/f {x,y', cb/f iz,t)) Ib/f iz,t), where the transverse ground- 
state Gaussians, nB/F,with widths ctb/f, and the axial functions, /b/f, are normalized 
to 1 and Nb/f respectively: 


vb/f {x,y;aB/F {z,t)) 


1 / -I- 1/^ \ 

tt^/Vb/f y 2al^.p{z,t)) 


( 8 ) 


I 12 

For both species, we define the axial density as Hid^b/f = |/b/f| • The relevance of 
the use of the Gaussian form in the confined directions for bosons and fermions was 
demonstrated in in Refs. [44] and [45], respectively, by showing that the respective 
factorization produces GSs close to their exact 3D counterparts. As demonstrated below, 
the Pauli principle makes the transverse width of the Fermi gas larger, in comparison 
with the Bose gas in the presence of a similar external potential. 

Replacing Tb/f by the factorized density in the Lagrangian density and integrating 
in the transverse plane of {x,y), the expression for action (1) is cast into takes the 
following form: 


5 


dtdz (Did,B + 'ClD,F + '^ID.Bf), 


where the effective ID densities are 

'^1D,B = (/B^t/n — /n^t/B) ~ ^1D,B’^1D,B — ClD.B, 


(9) 


( 10 ) 
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'^1D,F — (/p^t/p — /p<94/f) — ViD,pnid,p — eiD,p, 


( 11 ) 


- _ 1 5'bp 

'^1D,BP —-2~;-2 

TT Cr| + (Tf 


( 12 ) 


and eiD,B and eiD,p are the ID energy densities of the boson and fermion species, 
respectively: 


eiD,B 


2m-Q 


I^z/bI^ + 


9b h? 

47rcr| 2mBcr| 




n-iD.B, 


(13) 


eiD,p — 




2\2m-p 

, se 


5'P 




n 5/3 

n'iD,pn;^F) F) 


2mp 54 ^^ 


(14) 


with CiD,p = (3/5)(67r(2sF + 1))^'^^ in Eq. (14). Next, varying action S given by 
Eq. (9) with respect to and /p yields the respective Enler-Lagrange eqnations, i.e., 
the motion equations for the BFM in the ID approximation; 


ihdtfB 




2mB 


dl 




1D,B 


5'BP 


TT CTr + (T^ 


■I/fIV 


9b 

27rcr| 


I/b 


2mB(T^ 


-msWiBCTB 


/b, 


(15) 


i^dtfp = 
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rj2 I T/ I ^ hBF |„ ,2 

-- dz + Eid,F H-——^ /b 

2A2mF TT (Tg + (Tf 


I/p/ 


27rcr| 


ID 


Cf, 

2 mF 4/3 


l/i 


14/3 


1 2 2 I- 


(16) 


Relationships between ctb/f and /b/f are obtained by varying the ID action (9) with 
respect to (Tb/f: 


4 9b . 

X/,B^B-—niD,B — 0, 

ttib 27r 

4 2 „ 2/3 2 /: 

X/.Ft^B ~ y—sEpponiF, f'^f 

0 TTl'p 




\2mB 271 


9f _ „ 

—IJlD.F - 0, 


(17) 

(18) 


where X;,B/F = ™B/pWty/F “ 2gBFiiD,p/B/['F(<rB + <Ff)^1- 

Thus, Eqs. (15)-( 18) constitute a system of four ID coupled equations produced by 
the reduction of the underlying 3D system (5), (6). A numerical solution of Eqs. (17) 
and (18) by means of the Newton’s method produces the dependences of the transverse 
widths (Jb and ctf on the ID densities, Uid.b and Uid.Fi which are shown, respectively, 
in the left and right columns of Fig. 1, for the BFM with the attractive and repulsive 
interactions (the top and bottom rows, respectively). In all the cases, the Fermi species 
has, generally, a greater transverse width than its Bose counterpart. 
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Figure 1. (Color online) Color-coded charts for the transverse widths of the bosonic 
and fermionic species, ctb and crp (measured in units of ym), as functions of the bosonic 
and fermionic ID densities, niD,B/F- (a) cb and (b) ap for obf = —25 nm. (c) ctb and 
(d) (Tp for obf = 25 nm. The other parameters are ob/f = 5 nm and uJt, b/f = 1000 
Hz. Only fully polarized fermions are considered here. 


We consider for the bosonic density a range greater than for its fermionic 
counterpart because our calculations correspond to BFMs in which the number of bosons 
exceeds the number of fermions. For the attractive mixture (obf = —25 nm), it is 
observed that, for the bosons (see Fig. 1(a)) the width, ctb, slightly increases with niD,B 
and niD,F- For the fermions (see Fig. 1(b)), the width, uf, strongly increases with 
niD,F 5 due to the Pauli principle, and decreases with niD,B, because of the attractive 
boson-fermion interaction. In the repulsive mixture (obf = 25 nm), the situation for 
bosons is almost the same as in the case of the attraction, with the difference that ctb 
slightly decreases with the fermionic density, niD,F- For the fermions (see Fig. 1(d)) 
the increase in the boson density (niD,B) generates a signihcant increase in ap when 
the fermionic density (niD,F) is low. These four diagrams demonstrate the signihcant 
dependence of the widths of both species on the densities, thus showing the importance 
of treating ctb and ap as the variational variables. 

The above results can be explained by the dependence of the interaction between 
the species on their densities. To do this, we dehne a single scattering parameter. 


1 ^1D,F 




9e,B — QB 


(19) 
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and its counterpart for the fermions, obtained by replacing subscript B with F. Making 
use of the so dehned interaction parameter 19, in Eqs. 15 and 16, these equations assume 
the form given by 


and 


ihdtfB = 


i^dtfp = 




2mi 


d% + lhD,B + 


^?e,B I ^ 12 I 




2ms ah 








2 X 2 X 0 ^ 
C'f.id 


9% + ^1 


2mF (j. 


4/3 
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l/l 


1D,F 


14/3 


+ 


l?e,F I ^ |2 I 




2m-FX2(Tp 


-rripUJ-j-pO'-p 


/f- 


( 21 ) 


where the Eqs. 20 and 21 represented the motion equations for bosons and fermions 
separately, with effective scattering parameter ^fe^B and ge,F respectively. Figure 2 shows 
the variation of the effective scattering parameters, and ge^p, with the inter-species 
scattering length, obf- Figure 2(a) shows that the interaction between the bosons is 
suppressed when the interaction is attractive, and increases when it is repulsive. When 
the fermionic density is much lower that the bosonic density (circles), the influence of 
the variation of the interaction between the species on the effective interaction strength 
is minimal. Figure 2(b) shows the variation of the effective scattering parameter of the 
fermions {ge,F) under the same conditions as in Fig. 2(a), demonstrating the strong 
influence of the inter-species interaction on the fermions, due to the higher bosonic 
density. It is worthy to note that, for the range of obf presented in Fig. 2(b), ge,p is 
about three orders of magnitude higher than gp. The inset in Fig. 2(b) additionally 
shows the dependence of ge,p on the fermionic density niD,F (keeping niD,B fixed). A 
much larger number of bosons than that of fermions (the situation that that we address 
in this work), in the case of an attractive mixture, implies that the fermionic profile will 
be close to that of the bosonic gas. Accordingly, the Gaussian ansatz is appropriate for 
in both species, which is confirmed below by the comparison between 3D numerical and 
variational solutions. 


2.2. The two-dimensional variational system 

To derive 2D equations for the disc-shaped gas, the shape of the confinement potential 
is taken as (cf. Eq. (7) for the cigar-shaped configuration): 

Vb/f (r) = D 2 d,b/f {x, y) + -'mB/F^^l -QipZ^, (22) 

where the second term corresponds to the strong harmonic trap acting along 
the 2 ; direction. The corresponding factorized ansatz is adopted as Tb/f = 
4>b/f {.x,y)uB/F {z-,^b/f I where ub/f and 0 b/f are normalized to 1 and A^b/f, 

respectively, and ub/f is represented by the Gaussian wave function, 

.v2 


Ub/f {z; Cb/f) — — 

TT 
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1/4^1/2 


■exp 


B/F 


2^2 


(23) 


B/F 















Correlations and synchronization in a Bose-Fermi mixture 


9 




Figure 2. (Color online) The variation of the effective scattering parameters: (a) 
ge,B/gB and (b) ge.p/gp as a function of the interspecies scattering length ubf/cib for 
three values of niD.F/?^iD,B- The parameters are Wt b/f = 1000 Hz, ub/f = 5 nm and 
?T-iD,B = 100. The inset of the figure corresponds to ^bf/sb/f = 1 for two values of 

?T-1D,B- 


^B/F {,x,y,t) being the widths of the gas in the conhned direction. Substituting the 
factorized ansatz (23) into action (1) and integrating over z, we arrive at the following 
expression for the effective 2D action: 


<S — dtdxdy {£- 20 ,b + £ 2 F>,f + '^ 12 D,bf), 


h 

£ 2 F>,B = i— { 4 'B^t 4 >BB — “ h2D,B’^2D,B “ 620,B, 

£-21),F = * 77 ^ — ((>F<9t0p) — V 2 D,F’^ 2 D,F “ e 2 D,F, 

Za\ 


£ 


9bf 


2D,BF -^1/2 _ /p2 _L f2 


:? 2 ' 2 D,B? 2 . 2 D,F, 


(24) 

(25) 

(26) 
(27) 


'n/'Cb + 'Cf 

where n 2 D,B/F = |0b/f are the 2D particle densities of the boson and fermion 

species, and e 2 D,B and 620 ,f are their energy densities: 

Qb 1 9 / N 

+ 7 ''^bi^z,b'Cb ’ 2 - 20 ,b, (28) 


IV7 |2 I 

620,B — 77- V_L0B + 

zm-Q 


.\/^6 


-«2D,B + c2 
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620,F 




2\2m-p 


|Vb 0 f|' + 


9f 


2mT. 




n 2/3 
C2D,F'^2D F 


^ 2 D,F + 




4A2mF^F 4 




n2D,F, 


(29) 


with C'ao.F = (3/5)V2(6/(2 sf + l))'/=' 7 r in Eq. (29). 

The motion equations of the 2D system are obtained by the variation of the action 
S given by Eq. (24) with respect to variables 0 b and 0 f: 
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ihdtcpB = 




2m B 


72 I T/ I ^ ^bf 


F^2D,F 


9b 




■| 0 I 




(30) 


i^dt(t)F = 
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5'f 


2A2mF vr / \/^^f 
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2 m 




2 D,f| 0 f|^^^ + 




4A2mF^F 


+ -mFCU^ fCf 


0F. (31) 


Relations between ^b/f and 0b/f are prodnced by the Euler-Lagrange equations 
associated to ^b/f: 




^/,F'^F “ 


9b 


?^ 2 D,B^B 




mB 


= 0 , 


2^2 

Smp 




9f 


■■'n2D,F^F — 




= 0 


(32) 


(33) 


\/2^ A2mF 

where k/,f = mFCU^F + 2fi'BFn2D,B/[7r^'^^(.^| + 'Cf)^^^]- Equations (32) and (33) for ,^b/f 
were solved numerically. The four graphs in Fig. 3 show the dependence of .^b and .^f 
( in the left and right columns, respectively) on n 2 D,B/F 5 for the attractive and repulsive 
mixtures (top and bottom panels, respectively). The analysis for the attractive case 
can be divided into two parts: when the boson density is high, the fermions gas tends 
to compress, so that at low densities its width is found to be even less than that of 
the boson gas, a situation that tends to reverse for high fermionic densities; for low 
bosonic densities, the width of the boson gas rapidly decreases with the increase of the 
fermionic density, until becoming nearly constant, while, on the contrary, the width of 
the fermion gas increases progressively when its density is higher. Conversely, when the 
mixture is repulsive, the situation is similar to that outlined above for the ID setting: 
low densities of the mixtures cause the gas to compress, while high densities cause the 
mixture to expand, this trend being much stronger for the fermion component. 

To determine the influence that one species exerts on the other in 2D, we dehne 
the effective scattering parameters, cf. a similar dehnition (19) adopted in the ID case: 

2 \ “ 1 / 2 ' 

1 + (1 + dL) ] , (34) 

9b n2F>,B \ V'Cb, 


9e,B — 9B 


ge,F being produced by replacing subscript B with F. As shown in detail below, our 
calculations were performed assuming the same order of magnitude for the number of 
particles for both ID and 2D cases (for this reason that the effective ID density is 
about 10 times the magnitude of its 2D counterpart). With regard to this, conclusions 
concerning the effective interaction strength, which follow in the ID setting from Fig. 
2, apply to 2D case as well. 
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Figure 3. (Color online) Color-coded charts for the transverse widths of the bosonic 
and fermionic species, and (measured in units of ym), as functions of the bosonic 
and fermionic 2D densities, n 2 D,B/F- (a) and (b) ^p for obp = —25 nm. (c) ^b and 
(d) ^p for obp = 25 nm. The other parameters are ob/f = 5 nm and w^^b/f = 1000 
Hz. Only fully polarized fermions are considered here. 


3. Numerical Results 


Our analysis was developed for the GS and dynamics of perturbations around it in the 
presence of the 3D harmonic-oscillator trap of the following form: 

^b/f (r) = ^rrtB/F + ^^,b/f^^) , (35) 

In particular, for the GS of the ID and 2D settings we focused on determining the 
spatial correlation Cg between the spatial particle densities in both species, dehned as 


Cs {nB,np) 


(nBUp) 

V(nl) {nl) ’ 


(36) 


where u-b/f = np/p — (tt-b/f), () standing for the spatial average. For dynamical 
perturbations around the GS, a spatiotemporal correlation, which is dehned by replacing 
the spatial average with the spatiotemporal average, is known as the Pearson coefficient 
Cg-t [54]. The GSs in the ID and 2D cases were found by means of the imaginary-time 
integrations of Eqs. (15), (16) and (30), (31), respectively. 


3.1. Results in one dimension. 






















[10^ particles/m] 
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Figure 4. (Color online) Profiles of the particle density (as produced by the VA and 
3D calculations, see inset in panel (b)), and the width in the confined direction, for the 
bosonic and fermionic components, for five values of obf, as defined in the panels, (a) 
(b) uiD.F, (c) ctb, and (d) ctf- The parameters are Nb = 5x 10"^, A^f = 2.5x 10^, 
“b/f = 5 nm, b/f = 30 Hz and Wj b/f = 1000 Hz. The inset in (a) shows the 
difference between both methods, VA and 3D, by means of Auid.b = Pid,b — "raiD.B- 



Figure 5. (Color online) One-dimension overall percentage error for both species (see 
the definition in the text). The parameters are the same as in Fig. 4 
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3.1.1. Accuracy of the variational method as a function of scattering parameter obf for 
the ground states. As a starting point in the ID case with free coordinate z, we analyze 
the effect of the magnitnde and sign of the interaction parameter on the spatial prohle 
of both species, and the accnracy of the variational method (see Eqs. 15 - 18) compared 
to the 3D solntion, by varying the scattering length aD,B- The plots in Fig. 4 shows this 
sitnation, with panels 4(a) and 4(b) corresponding to the prohles of niD,B and niD,F, 
respectively, where the inset in panel 4(b) specihes whether the prohle was prodnced by 
the VA or by means of the 3D solntion, while 4(c) and 4(d) correspond to the prohles 
of (Tb and ap. The mixtnre with many more bosons than fermions is considered here: 
A^b = 5 X 10^, A^f = 2.5 X 10^. Because of this condition, the bosonic prohle is mainly 
determined by its self-interaction and the external potential. First, we deal with the 
inhuence of the interspecies interaction (^^bf) on the density prohle, and then address 
an error resulting from the variational approach. Figure 4(a) shows that variations of 
the bosonic density prohle are very small in comparison to the signihcant changes of the 
inter-species scattering length. The situation is opposite for the fermionic species. As the 
repulsive scattering length increases, the fermions tend to be pushed to the periphery of 
the bosonic-gas density prohle. This phenomenon is known as demixing [27, 32, 35, 52]. 
On the other hand, for the attractive case, fermions are, naturally, concentrated in 
the same region where the bosons are located. Figures Fig. 4(c) and 4(d) show that 
the width of the bosonic prohle (Fig. 4(c)) slightly increases as going from the inter¬ 
species attraction to repulsion. A similar trend is observed for fermions, see Fig. 4(d), 
but amplihed in the spatial zone of the interaction with the bosons, where the gas is 
compressed in the case of the attraction and expands in the case of the repulsion. It is 
noted that the fermionic component expands in the conhned direction much more than 
its bosonic counterpart, and that the fermionic width huctuates markedly with changes 
in density. Now, to analyze the accuracy of the method we solved Fqs. 5 and 6 and 
obtained the ID density associated through of Pid,b/f = f f |Tb/f| dxdy. The inset 
in panel 4(a) shows that the difference between the bosonic prohles obtained by both 
methods (VA or 3D) is ~ 2% of the maximum density for all cases (the fact that the 
error changes very little with variations in obf is a consequence of the greater number of 
bosons). Figure 4(b) shows that, for the case of the attractive mixture, the variational 
prohle is very close to the 3D result, in particular for the case of obf = —Qnm. For the 
repulsive mixture, it is observed that the error increases, which is a consequence of the 
lower fermionic density at the center of the 3D harmonic potential, dominated for the 
bosons, hence a monotonously decreasing function in the transverse direction, such as 
the Gaussian, is not a sufficiently good approximation. We dehne the global error for the 
variational method as E%^ib, = f |piD — ''^idI dz (for both species). Figure 5 shows the 
global error for both species as a function of the scattering parameter, obf- This picture 
demonstrates that the error for the bosonic species is around 2%, and it does not change 
much (as already noted in the comment to Fig. 4(a)). For the fermions, the error attains 
a minimum value ~ 0.5% at obf ~ 4nm, thus corroborating that the Gaussian is a very 
good approximation. This is a consequence of the fact that, for this value of obf, the 
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Figure 6. (Color online) The color-coded plot of the spatial correlation Cg, defined 
as per Eq. (36), as a function of the fermionic number, N-p, and the scattering length 
of the inter-species interaction, obf- Other parameters are N-q = 5 x 10^, ob/f = 5 
nm, LOz.B/F = 30 Hz and Wt b/f = 1000 Hz. 


interspecies interaction term practically compensates the Pauli repulsion term, making 
the dynamics of the fully polarized Fermi gas close to that giverned by the Schodinger 
equation (recall that the Gaussian is the solution for the ground state). When the 
mixture becomes more attractive, the fermionic dynamics is dominated by the bosons, 
producing a similar error for both species, while for the repulsive mixture the Gaussian 
approximation is not appropriate. Interestingly, for the non-interacting mixture, the 
error for the fermions is smaller than for the bosons, because the fermionic density is 
very low, making the self-interaction terms weak in comparison to the external potential, 
therefore it is appropriate to use the Gaussian ansatz for describing the ID dynamics. 
Using the TF prohle may produce, in principle, better results, but the respective analysis 
would be cumbersome. 

3.1.2. Spatial correlations in the ground state. As the hrst application, we calculate the 
spatial correlation between densities of the two species (G^), as dehned by Eq. (36) for 
a wide range of parameters, ranging from very attractive to to very repulsive mixtures. 
We keep in mind that the variational calculation produces a minor error for attractive 
mixtures, therefore the values obtained for the repulsive case will be considered as help to 
asses the validity of our analysis. The eventual objective is to produce a parameter which 
determines the mutual influence of both species. In the Fig. 6, the spatial correlation Cg 
is shown versus N-p and obf , which shows that, for an attractive interaction, Cg reaches 
values greater than 0.8 (in the mixed state), whereas for the repulsive interaction Cg 
decreases until reaching values around Cg = 0.3 (in the demixed state). It is also 
noted that the number of fermions strongly affects the dependence of the correlation 
on the inter-species scattering length. When Np increases, the difference between the 
maximum and the minimum value of Cg decreases, as obf varies from —25 nm to -|-25 
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Figure 7. (Color online) Space-time diagrams for the conservative dynamics of bosons 
(above) and fermions (bottom), for three different values of the interspecies scattering 
parameter: (a)—(b) obf = —6nm, (c) —(d) obf = —12nm, and (e) —(f) obf = —llnm. 
The initial conditions are the same in all cases and are described in the text. Other 
parameters are the same as in Fig. 4 

nm, tending to reach valnes close to Cs = 0.8, which is the valne attained in the absence 
of the interaction, depending only on the presence of the same trap acting on both 
components. 

3.1.3. Dynamics near the ground state We begin the stndy of conservative dynamics, 
considering a mixture with arbitrary initial conditions for the ID helds: we assume 
a Gaussian shape along the z-a.xes with widths (standard deviation) of lOO/im and 
200/im for bosons and fermions, respectively. Figure 7 shows three cases of the 
temporal evolution with these initial conditions for obf = —Qnm, obf = — 12nm, 
and Obf = — Idnm. In the hrst case (panels 7(a) and 7(b)), it is observed that 
the densities converge towards the center of the potential in an expected pattern of 
oscillations around the potential minimum; in addition, the fermions are affected by 
bosons, as can be seen for the mark left by the bosons in the fermionic density. For the 
second case (panels 7(c) and 7(d)), it is observed that the increase in the magnitude 
of the attractive interaction generates dark solitons in the fermionic density, some of 
which shows oscillatory dynamics very similar to that observed, experimentally and 
theoretically, in Refs. [55]-[60]. Finally, the last case (panels 7(e) and 7(f)) shows 
that further increase in the magnitude of the interspecies interaction generates a larger 
number of dark solitons. In other words, we show that the attractive interaction of 
fermions with bosons in a state different from the ground state eventually generates a 
gas of dark solitons. 
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Figure 8. (Color online) (a) and (b) Spatiotemporal evolution of the bosonic and 
fermionic densities, respectively, in the case of suddenly switching the inter-species 
interaction in the ID ground state which corresponded to obf = 0. (c) and (d) Bosonic 
and fermionic density profiles at the end of the simulations (the blue lines) compared 
to the ground state (red line), at the same values of physical parameters, which are 
Nb = 5 X 10"^, N-p = 2.5 X 10^, Ob/f = 5 nm, obf = —15 nm, b/f = 30 Hz, and 
iOi B/F — 1000 Hz. 


Now we aim to address the system’s dynamics in the vicinity of the GS. We start 
with the GS found in the absence of the inter-species interaction (obf = 0). Then, at 
t = 0, we switch the interaction on, which may imply the application of the magnetic 
held, that gives rise to obf 7^ 0 via the FR. Figure 8 shows an example of the ensuing 
evolution in the attractive mixture with obf = —15 nm, the fermionic component being 
fully polarized. Panels 8(a) and (b) display space-time diagrams of oscillations of the 
bosonic and fermionic densities in the dynamical state. Panels 8(c) and (d) compare the 
density prohles at the end of the simulations with the GS of the interacting mixture. It 
is seen that the dynamics amount to rapid huctuations around the GS. Figure 7 makes 
it clear that the spatiotemporal dynamics of fermions in panel 8(b) does not reduce to 
noise around the ground state. It rather corresponds to a gas of dark solitons oscillating 
in the harmonic potential. Moreover, the fact that the spatiotemporal diagram has 
been produces for the time when the system is in the stationary regime suggests that 
the dark solitons are fully stable. Note that the noise observed in the bosonic component 
is clearly caused by the fermionic dark solitons. 

Now, we measure how correlated the densities remain in time. Similar to what is 
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Figure 9. (Color online) The evolution of spatial correlation Cs in the dynamical 
states near the ground state for five values of obf, as indicated in the figure. The 
parameters are the same as in Fig. 8. 


observed in Fig. 6, Fig. 9 shows that the spatial-correlation parameter Cg, dehned as 
per given Eq. (36) rapidly decreases following the transition from negative to positive 
valnes of obf (the mixing-demixing transition). When the inter-species interaction is 
more attractive, it increases the correlations, imposing spatiotemporal synchronization 
of density flnctnations. In the case of strong repulsion, obf = 50 nm, correlation 
fluctuations are suppressed because the system falls into a demixed state. 

Figure 10 presents the analysis of the spatiotemporal synchronization of the 
mixture, where three possible regimes are considered for the fermions: fully polarized, 
BCS, and unitarity. Parameters of the the Lagrangian density for each fermionic regime 
are given in Table 1. The synchronization was calculated through the Person coefficient 
Cg-t dehned above. For each of the three cases, the spatial correlation, Cs, was 
calculated too. When the interaction is attractive, there is no discrepancy between 
the synchronization and correlation curves. On the other hand, in all the three regimes 
it is observed that, as the interaction becomes more repulsive, the synchronization of 
the mixture decreases in comparison to the spatial correlation. This trend enhances as 
the system approaches the demixing transition. Another generic feature is that stronger 
correlations between the species near the GS imply a stronger synchronization on the 
spatiotemporal dynamics too. 

The correlation properties are different for the distinct fermionic regimes. Cg^t is 
always positive for the system with fully polarized fermions and in the BCS regime, 
decreasing as the interaction gets more repulsive. On the other hand, the unitarity 
regime features signihcant differences: at hrst, Cg-t grows to a maximum value at 
obf = —10 nm; then it decreases, reaching negative values of Cg^tjCg in the demixed 
state. In all the cases, no signihcant difference are observed between Cgt and Cg in the 
case of the attractive Bose-Fermi interaction, i.e., a highly correlated GS supports the 
dynamical synchronization too. 
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Figure 10. (Color online) Spatiotemporal correlation Cs-t versus obf for the 
ground state (GS) and in the dynamics near the ground state (“DNGS”) for three 
fermions regimes: polarized, BCS, and unitarity. The parameters are A^b = 5 x 10^, 
Np = 2.5 X 10^, cib/f = 5 nm, uJz^B/F = 30 Hz and uit^B/F = 1000 Hz. Table 1 provides 
parameters that correspond to the different regimens in the Lagrangian density. 


3.1.4- Accuracy of the variational method for the spatiotemporal dynamics In the 
Sec. 3.1.3 some examples of spatio-temporal dynamics were presented without discussing 
the accuracy of the results. In Fig. 11 the spatio-temporal dynamics of the ID density 
is displayed as produced by the solution of the 3D equations 5 and 6, using the same 
procedure as in Section 3, for initial conditions similar to those used in Fig. 7 but with 
Obf = — lOnm. The initial conditions for the 3D dynamics are given by the ansatz in 
Eq. 8, with the Gaussian prohle along the z-axes, similar to those used in the Fig. 7. 
Panels 11(a) and 11(b) show spatio-temporal diagrams of the bosonic and fermionic 
density, respectively, making the emergence of dark solitons obvious, see also Fig. 7. 
This result corroborates that the dark solitons will emerge too in the 3D dynamics, 
which is approximated by the present ID model. The other panels show a comparison 
of the ID spatial prohles, as obtained from 3D simulations and the ID VA, for three 
instants of time: t = 0ms (panels 11(c) and 11(d)), t = 25ms (panels 11(e) and 11(f)), 
and t = 50ms (panels 11(g) and 11(h)). The results demonstrate that the VA prohles 
are very similar to their counterparts produced by the 3D simulations, hence the present 
approximation provides good accuracy and allows one to study dynamical features of 
the Bose-Fermi mixture in a sufficiently simple form. 

3.2. Results in two dimensions. 

3.2.1. Accuracy of the variational method as a function of scattering parameter obf for 
the ground state. In the 2D setting, we studied the GS of the BFM, hxing the number 
of bosons to Nb = 5 x 10“^ (similar to the value used in the ID case), which was and 
much greater than the number of fermions. Fig. 12 shows the prohle of the 2D bosonic 
and fermionic densities, n 2 D,B/F (top), as obtained from the VA and 3D simulations, and 
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Figure 11. (Color online) Comparison of the dynamics, as obtained from the ID 
variational approximation, and from the 3D simulations. Spatiotemporal diagrams for 
bosons (a) and fermions (b) are obtained from the 3D simulations. The other panels 
show spatial profiles for: (c) —(d) t = 0ms, (e)—(f) t = 25ms, and (e)—(f) t = 50ms. 
Here obf = — lOnm, the initial conditions and other fixed parameters being the same 
as in Fig. 7. 


the widths of the mixture in the transverse direction, ^b/f (bottom), with respect to the 
radial coordinate. To obtain the 2D prohle from the 3D simulations (results of solving 
Eqs. 5 and 6), we integrate the 3D density along the z so that P2D,b/f = / |^2D,b/f| dz. 
The panels for the bosonic and fermionic components are displayed on the left and right, 
respectively. The results are very similar to what was observed in the ID case in Fig. 4, 
for the same value of the physical parameters. Therefore, the following conclusions 
are also similar to what inferred in the ID setting: the repulsive mixture concentrates 
the bosons at the center, while the attractive mixture concentrates both species at 
the center. Panels 12(c) and 12(d) show that only the width of the fermionic density 
prohle varies signihcantly with the change of the scattering length of the inter-species 
interaction, which is consequence of a greater number of bosons than fermions. It is 
clearly seen that fermions are stronger conhned when the interaction is attractive, and 
their spatial distribution signihcantly expands when the interaction is repulsive. Similar 
results have been observed in [27, 28, 32]. 

Now, to compare the results obtained from the VA with those produced by the 3D 
simulations, we note that both prohles are practically identical, except for the repulsive 
case in which a discrepancy is observed (smaller than in the ID case). The inset in 
panel 12(a) shows that the diherence between the two results has a magnitude of nearly 
three orders of magnitude lower than the density itself. The error is lower than in the 
ID case because the 2D reduction is closer to the full 3D model. 

Now, similar to the ID case, we dehne the overall percentage error of the VA as 
A'%, 2 D = / / |P 2 D — ''^ 2 d| dxdy (for both species). Figure 13 shows the error for both 
species as a function of interspecies scattering parameter (obf)- For bosons it is around 
of 0.2% (one order of magnitude lower than in the ID case), and does not change much 
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Figure 12. (Color online) The radial profile of the 2D particle density (for the VA and 
3D simulations, see the inset in panel (b)), and the respective width in the transverse 
direction (z) of the bosonic and fermionic profiles, for five values of obf (see insets in 
the panels): (a) n 2 D,B, (b) n 2 D,F, (c) and (d) .^f- The parameters are Nq = 5 x 10"^, 
TVf = 2.5 X 10^, ob/f = 5 nm, iOz , B/F = 1000 Hz, and uJx,b/f = w^^b/f = 30 Hz. The 
inset in panel (a) shows the difference between both methods, VA and 3D, where 

An2D = P2D — n2F). 


(as shown in the inset to Fig. 12(a)). For fermions the error is greater than for bosons 
thronghout the observed range, bnt it is qnite small for the attractive mixtnre. Like 
in the ID case, the error increases for the repnlsive mixtnre, bnt remains lower than 
that calculated for bosons in the ID case. In summary, the 2D approximation is very 
accurate, implying, as shown below, that the description of the dynamics is also more 
accurate in this case. 

3.2.2. Spatial correlations of the two-dimensional density Here we again apply the 
formalism to the calculation of the spatial correlation Cg. The Fig. 14 shows the 
dependence of the spatial correlation C*, defined as per Eq. (36), on the scattering 
length obf for three values of N-p and fixed Np = 5 x 10^, in the regime of the fully 
polarized fermionic component. The case with a smaller number of fermions {N = 2500, 
the same as in Fig. 12) gives rise to the greatest contrast in the values of the correlation; 
from a value close to 1 for g-Qp = —25 nm to a value near zero for g-pp = 25 nm. This 
transition is more drastic than that observed for the same parameters in the ID case. 
This is mainly due to different factors multiplying the interaction terms in the 2D 
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Figure 13. (Color online) 2D overall percentage error for both species (the definition 
is given in the text). Parameters are the same as in Fig. 12. 



Figure 14. (Color online) Spatial correlation Cg of the ground state of the 2D mixture 
as a a function of obf, for three values of the number of fermions. The other parameters 
are: A^b = 5 x 10^, ob/f = 5 nm, w^j^b/f = ^y,B/F = 30 Hz and 0Jz,b/f = 1000 Hz. 
Fermions are in the fully polarized state. 


equations (31), ( 32), in comparison with their ID counterparts (15), (16). The contrast 
between the mixing and demixing decreases as the number of fermions increases, due to 
the strengthening of the self-interaction of the fermions, hence the correlation gets less 
sensitive to the inter-species scattering length, obf- 

The Fig. 15 shows the spatial correlation, dehned as per Eq. (36), for the three 
fermionic regimes: polarized, BCS and unitarity. The hgure is a counterpart of Fig. 10 
drawn for the ID case, but only for the spatial correlation in the GS {Cg)- One may 
conclude that these fermionic regimes show the behavior similar ti that observed in 
the ID case, a difference being that the three curves demonstrate stronger demixing 
when obf changes from positive to negative values. In the unitarity regime it is again 
observed that the correlation reaches a maximum close to 1 at obf ~ —10 nm, dropping 
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Figure 15. (Color online) Spatial correlation Cg of the ground state of the 2D 
mixture as a a function of obf, for three fermionic regimes: polarized, BCS, and 
unitarity. The fixed parameters are: Nb = 5 x 10^, iVp = 2.5 x 10^, ob/f = 5 nm, 
<j-’x,B/F = Wy,B/F = 30 Hz and w^^b/f = 1000 Hz. 


to negative values when the mixture is strongly repulsive. 

3.2.3. Accuracy of the variational method for the spatiotemporal dynamic Figure 16 
shows, for the attractive mixture (^^bf = — 12nm), the comparison between the results 
obtained by the t3D integration (Eqs. 5 and 6) with the results obtained by the 2D 
integration of the variational equations (Eqs. 30 and 31). The approach is similar to 
that developed in Sec. 3.1.4 for the ID case. It starts from arbitrary initial conditions 
in 2D, viz., 2S fields with Gaussian shapes in the radial direction, a standard deviation 
of 50/im for bosons and lOOpm for fermions being assumed. The initial conditions 
for the 3D dynamics are obtained with the ansatz given by Eq. 23 and assuming a 
Gaussian in the x|/-plane. Panels 16(a) and 16(b) display two spatiotemporal diagrams, 
obtained from the 3D integration, showing a cut along the x axis of the 2D density; since 
the dynamics, for the time range observed, obeys the cylindrical symmetry, this cut is 
representative. It is seen that the boson profile clearly affects the fermion prohle. As 
the initial conditions do not correspond to the ground state, the fermions are attracted 
to the center by the bosonic pattern and are subject to the conhnement imposed by the 
harmonic potential, generating strong interaction with the boson core. A gray soliton 
ring is generated (in panel 16(b) it is seen as a gray soliton). In panels 16(c) and 16(d), 
the radial profile obtained by the 3D integration and by VA for three instants of time 
are compared, showing that it is not possible to distinguish them. The inset in panel 
16(c) shows that, in the observed range, the 2D overall error is less than 1% (at initial 
time, it is zero because the 3D initial conditions are obtained from the 2D conhguration 
using the ansatz, Eq. 23). The approach delivers very good results. In addition, in 
panel 16(d) it was observed that the rings are gray (the minimum value of the fermion 
density is not zero). Panels 16(e) and 16(f) show diagrams of the 2D density in the 
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Figure 16. (Color online) 2D dynamics, as produced by the 2D VA, and by the 
3D simulations for the attractive mixture of bosons (above) and fermions (below). 
Spatiotemporal diagrams (a) and (b), produced by the 3D simulations, show how the 
profiles of the densities of particles in the x direction evolve in time (it is a cut of 
the 2D density with cylindrical symmetry). The curves in panels (c) and (d) shows 
the radial profiles for three different instants of time (see the inset in panel (d)), as 
produced by both the VA and 3D simulations. Diagrams (e) and (f) show 2D densities 
in the xj/-plane for the last moments of time in panels (a) and (b). The inset in panel 
(c) shows the evolution of 2D error (A% 2 d) in time. Initial conditions are described 
in the text, and parameters are the same as in Fig. 12 with obf = —12nm. 


a;t/-plane for the last moment of time, highlighting the gray ring soliton observed in the 
fermionic component, which expands wherein the boson density is zero. We note that, 
in this case, the approximation is very good and allows one to considerably reduce the 
computing time needed to observe the 2D dynamics. The dynamics of the mixture thus 
produces a vast variety of localized nonlinear structures. 

4. Conclusions 

The incorporation of the width of the conhned species as a variational variable shows 
that the term accounting for the interaction between the bosonic and fermionic species 
in the lower-dimensional equations is inversely proportional to the width in the conhned 
direction(s). The derived systems of equations, which include the algebraic equations 
for the widths as functions of the ID or 2D bosonic and fermionic densities. These 
approximations yield solutions which agree well with the full 3D theory. It is worthy 
to note that the ensuing relation between the widths and densities signihcantly varies 
with the strength of the inter-species interaction, which can be tuned by means of the 
Feshbach resonance. For very low densities, our equations are similar to the simplest 
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models which postulate constant conhnement width in the transverse direction(s). As an 
application of the variational method, we have studied the spatial correlation between 
the species, in which a strong dependence on the interaction parameter was observed. 
This sensitivity becomes stronger for the repulsive interaction, which causes demixing 
of the species. 

It was observed too that, in ID and 2D alike, the spatial correlation between the 
species is also strongly affected by the inter-species interaction strength. In addition, 
it was concluded that the greatest contrasts in the correlation, changing from the large 
correlation in the highly attractive mixtures to a weak correlation in the strongly 
repulsive ones, occur in the unitarity fermionic regime. Further, in the ID setting 
the dynamics near the ground state demonstrates that the spatiotemporal correlation 
follows the same trends as the spatial correlation, the largest (but still small) difference 
being observed for strongly repulsive mixtures 

Finally, it is relevant to mention that the mixture can also give rise to stable 
dark solitons in the Fermi field, in both the ID and in 2D cases. However, as the 
hydrodynamic conditions was assumed, the solitons may occur in some form but are 
unconfirmed and should not be believed in too strongly. A detailed characterization 
of these states and the validity of the hydrodynamic approximation for them will be 
presented elsewhere. 
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